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Abstract 

We derive algorithms to perform a spin spherical harmonic transform and inverse for functions of arbitrary spin 
number. These algorithms involve recasting the spin transform on the two-sphere as a Fourier transform on the 
two-torus T^. Fast Fourier transforms are then used to compute Fourier coefficients, which are related to spherical 
harmonic coefficients through a linear transform. By recasting the problem as a Fourier transform on the torus we 
appeal to the usual Shannon sampling theorem to develop spherical harmonic transforms that are theoretically exact 
for band-limited functions, thereby providing an alternative sampling theorem on the sphere. The computational 
complexity of our forward and inverse spin spherical harmonic transforms scale as 0{L^) for any arbitrary spin 
number, where L is the harmonic band-limit of the spin function on the sphere. The algorithms also apply for functions 
with arbitrary band-limit and not only powers of two. Numerical experiments are performed and unfortunately the 
forward transform is found to be unstable for band-limits above L :^ 32. The source of this instability is due to 
the poorly conditioned linear system relating Fourier and spherical harmonic coefficients. The inverse transform is 
expected to be stable, although it is not possible to verify this hypothesis. In ongoing work we are attempting to 
renormalise or reformulate the algorithms in such a way as to eliminate the numerical stability problem. 
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1. Introduction 



Data are defined inherently on the sphere in many practical applications, including computer graphics 
{e.g. [26]), planetary science {e.g. [32,37,38]), geophysics {e.g. [29,31,34]), quantum chemistry {e.g. [ , ]) 
and astrophysics {e.g. [3,4,14,15]), for example. In a number of these applications it is insightful to compute 
the spherical harmonic representation of the data. For example, observations of the anisotropics of the 
cosmic microwave background (CMB), that are inherently made on the celestial sphere, contain a wealth 
of information about models of the early Universe. The anisotropics of the CMB are seeded by primordial 
perturbations, which, in linear perturbation theory, evolve independently for each Fourier mode. Solutions of 
models of the early Universe are therefore solved more conveniently in Fourier space and manifest themselves 
in the shape of the angular power spectrum of CMB anisotropics. Cosmologists compute the angular power 
spectrum of observations of the CMB, through a spherical harmonic transform, to make contact with theory 
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and to estimate parameters of cosmological models {e.g. [10]). Recent and upcoming observations of CMB 
anisotropics made on the celestial sphere are of considerable size, with maps containing approximately three 
[15] and fifty [24] million pixels respectively. Consequently, the ability to perform a fast spherical harmonic 
transform for data of considerable size is essential in assisting cosmologists to understand our Universe. 

A number of works have addressed the problem of computing a spherical harmonic transform efficiently. 
We first survey the related literature for computing the scalar transform on the sphere, before discussing 
approaches for computing the harmonic transform of spin functions on the sphere. By performing a simple 
separation of variables the spherical harmonic transform may be rewritten as a Fourier transform and an 
associated Legendre transform. This already reduces the complexity of computation from 0{L^) to 0{L^)^ 
where L is the harmonic band-limit of the function considered. If this separation of variable is performed 
first, then to reduce the complexity of the spherical harmonic transform further reduces to the problem of 
performing a fast associated Legendre transform. 

First attempts to compute a scalar spherical harmonic through a fast Legendre transform were performed 
by Orszag [ ] and were based on a Wentzel-Kramers-Brillouin (WKB) approximation. Alternative ap- 
proaches using the fast multipole method (FMM) [ ] have been considered by Alpert & Rokhin [ ] and by 
Suda & Takani [ ]. All of these methods are necessarily approximate (for any choice of pixelisation of the 
sphere) and the complexity of these algorithms scales linearly with the desired accuracy. 

Other approaches based on the separation of variables have been developed for particular pixelisations 
of the sphere, such as HEALPix^ [12], IGLOO ^ [7] and GLESP^ [8]. In ah of these cases the resulting 
algorithms are of complexity 0{L^) and for HEALPix and IGLOO pixelisations only approximate quadrature 
rules exist. Due to the approximate quadrature the spherical harmonic transform algorithms of HEALPix and 
IGLOO are not theoretically exact; a forward transform followed by an inverse transform will not recover the 
original signal perfectly, even with exact precision arithmetic. Nevertheless, all of these pixelisation schemes 
and their associated harmonic transform algorithms are of sufficient accuracy for many practical purposes 
and have been of considerable use in the analysis of CMB data. 

A sampling theorem on the sphere was first developed by Driscoll & Healy [ ], providing a framework 
for a spherical harmonic transform that is theoretically exact. Moreover, Driscoll & Healy also present a 
divide and conquer approach to computing a fast associated Legendre transform in the cosine basis, using a 
fast Fourier transform (FFT) to project onto this basis. The resulting algorithm is exact in exact precision 
arithmetic and has computational complexity \0g2L)., but is known to suffer from stability problems. 

Potts [25] adapted the Driscoll & Healy approach to use a fast cosine transform (as suggested as future 
work in [9]) to reduce computation time (but not complexity) and introduced some stabilisation methods. 
Healy et al. [13] readdressed the work of Driscoll & Healy, reformulating the sampling theorem on the sphere 
and developing some variants of the original Driscoll & Healy algorithm. A fast cosine transform is used 
here and a number of alternative algorithms are presented, including the so-called semi-naive, simple split 
and hybrid algorithms. The semi-naive algorithm is a more stable variant of the original Driscoll & Healy 
algorithm, while the simple-split algorithm is more simple and stable but with an increased complexity of 
0{L^ log L). The hybrid algorithm uses both of these two algorithms and spits the problem between them. 
The hybrid algorithm appears to achieve a good compromise between stability and efficiency. However, the 
overall complexity of this algorithm is not clear since it depends on the spit between the semi-naive and 
simple split algorithms and on some user specified parameters. Implementations of these algorithms [ ] are 
available for download. Due to the divide and conquer approach first proposed by Driscoll & Healy for 
the fast associated Legendre transform, these algorithms [9,25,13] are all restricted to harmonic band-limits 
that are a power of two. 

An alternative algorithm is developed by Mohlenkamp [19,20], that is again based on a separation of 
variables and a fast associated Legendre transform. This method is based on the exact sampling theorem on 
the sphere of Driscoll & Healy and a compressed representation of the associated Legendre functions using 
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WKB frequency estimates. Rigorous bounds are derived for the number of coefficients required to achieve 
a given accuracy. The method is again necessarily approximate, but the approximation can be controhed 
independently of complexity. The complexity of this algorithm is 0{L'^ log 2L) and in practice very good 
accuracy and stability is achieved. However, the algorithm is also restricted to band-limits that are a power 
of two. 

Very little attention has been paid to fast algorithms to perform the spherical harmonic transform of spin 
functions. Spin functions on the sphere also arise in many applications. For example, the linear polarisation 
of the CMB anisotropics is described by the Stokes parameters Q and /7, where the quantity Q ^ ill is a 
spin ±2 function on the sphere [3n]. The reason little attention has been paid explicitly to the transform 
of spin functions may be because spin lowering and raising operators can be used to relate spin spherical 
harmonic transforms to the scalar transform. This is done explicitly and in a stable manner by Wiaux et al. 
[35] for the spin ±2 case. Scalar spherical harmonic algorithms developed by Healy et al. are then applied 
to yield a theoretically exact spin ±2 algorithm with the same complexity as the Healy et al. algorithm that 
is adopted. In general, this type of approach may be used to compute the spin transform for arbitrary spin 
(required recurrence relations are given in [00]), however the complexity of the resulting spin algorithm will 
be scaled by the spin number. 

In this work we derive a general algorithm to perform a spin spherical harmonic transform for functions 
of arbitrary spin number. The complexity of the algorithm is independent of spin number. Furthermore, any 
harmonic band-limit may be considered and not only powers of two. Our spin spherical harmonic transform 
algorithm essential involves recasting the spin transform on the sphere as a Fourier transform on the torus. 
FFTs are then used to compute Fourier coefficients, which are related to the spherical harmonic coefficients 
of the function considered through a linear transform. By recasting the problem as a Fourier transform on 
the torus we appeal to the usual Shannon sampling theorem to develop spherical harmonic transforms that 
are theoretically exact for band-limited functions, thereby providing an alternative sampling theorem on 
the sphere. The remainder of this paper is organised as follows. Mathematical preliminaries are presented 
in Sec. 2. In Sec. 3 we derive our spin transform algorithms and discuss a number of subtleties of the 
algorithms and practical considerations. In Sec. 4 we examine the memory and computation complexity of 
the algorithms and evaluate their stability with numerical experiments. Concluding remarks are made in 
Sec. 5. 

2. Mathematical preliminaries 

Before deriving our fast algorithms it is necessary to outline some mathematical preliminaries. Harmonic 
analysis of scalar functions on the two-sphere S^ and on the rotation group SO (3) is reviewed first, before 
discussing the analysis of spin functions on S^. By making all definitions explicit we hope to avoid any 
confusion over the conventions adopted. The reader familiar with harmonic analysis on S^ and SO (3) may 
choose to skip this section and refer back to it as required. However, in the final part of Sec. 2.2 we present 
some important symmetry relations of Wigner functions that are required in the derivation of our algorithms. 

2.1. Scalar spherical harmonics 

The spherical harmonic functions Y£m{^) with natural ^ G N and integer m G Z, |m| < ^ form an 
orthogonal basis in the space of L^(S^, d^}{u)) scalar functions on the two-sphere S^, where uo = (6>, (p) G S^ 
are the spherical coordinates with co-latitude G [0, tt] and longitude cp G [0, 27r) and dQ{u;) = sinOdOdcp 
is the usual rotation invariant measure on the sphere. The spherical harmonics are defined by 
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are the associated Legendre functions. We adopt the Condon-Shortley phase convention, with the (— 1)^ 
phase factor included in the definition of the associated Legendre functions, to ensure that the conjugate 
symmetry relation ¥^^{00) = (— l)^y£^_^(a;) holds, where the superscript * denotes complex conjugation. 
The orthogonality and completeness relations for the spherical harmonics defined in this manner read 

/ dn{u) Yim{(^) Yl;^;{ij) = 5it5mm' (3) 

and 

oo i 

^ YeM Yl^{uj') = S{co - u;') (4) 

£=0 m^-£ 

respectively, where 6ij is the Kronecker delta, 6{uj — uo') = 8{ip — ^p') 6{cos6 — cos 6') and 6{x) is the Dirac 
delta function. Since the spherical harmonic functions form a complete, orthogonal basis on the sphere, any 
square integrable function on the sphere / G L^(S^, drt{u)) may be represented by the spherical harmonic 
expansion 

oo i 

/M = E E fe^YeM, (5) 

£^0 m^-£ 

where the spherical harmonic coefficients are given by the usual projection on to the spherical harmonic 
basis functions 

fem= I dn{u;)f{u;)Y;M' (6) 

The conjugate symmetry relation of the spherical harmonic coefficients of a real function is given by 
fi^-m = (~l)^//m' which follows directly from the conjugate symmetry of the spherical harmonic func- 
tions. 

2.2. Wigner functions 

We describe here representations of the three-dimensional rotation group SO (3). Any rotation p G SO (3) 
is uniquely defined by the Euler angles p = (a, 7), where a G [0, 27r), /3 G [0, tt] and 7 G [0, 27r). We adopt 
the zyz Euler convention corresponding to the rotation of a physical body in a fixed coordinate system about 
the z, y and z axes by 7, P and a respectively. The Wigner i^-functions Df^^{p) are the matrix elements 
of the irreducible unitary representation of the rotation group S0(3). Moreover, by the Peter- Weyl theorem 
on compact groups the matrix elements Df^^ also form an orthogonal basis on L^(S0(3), d^(p)), where 
dg{p) = sin /3 da d/3 dj is the invariant measure on the rotation group. The orthogonality and completeness 
relations for the Wigner D-functions read 



f I Stt 

/ D^^^^p) Dil^\p) = - du'^mm'^nn' (7) 
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and 



oo . I I 

E ^ E E D'^nip) Dtnip') = 5{P - P') (8) 

£=Q m^-£n^-£ 

respectively, where 5{p — p') = S{a — a') 6 {cos P — cosP^)S{j — Y). Since the Wigner /^-functions form 
a complete, orthogonal basis on the rotation group SO (3), any square integrable function on the rotation 
group F G L^(S0(3), dg{p)) may be represented by the Wigner /^-function expansion 

oo £ £ 

P^P) = E ^ E P'ran DUp) , (9) 

£^0 m^-£n^-£ 
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where the Wigner harmonic coefficients are again given in the usual manner by 



F^n = / Mp) F{P) D'Z^ip) ■ (10) 

^SO(3) 

The Wigner /^-functions may be decomposed in terms of the reduced Wigner d-functions by 

Dl„{a,P,j) = e-™« dl^iP) e-»^ . (11) 
The real (i-functions may be defined using the relation 



dLniP) =(-l)«- + m)\{£ - m)\{£ + n)\{£ - n)\ 

/ ^\2£—m—n—2k / _\ m+n+2/c 

^ k\{£-m-k)\{£-n-t)\{m^n^k)\ ' ^ ^ 

where the sum is performed over all values of k such that the arguments of the factorials are non- negative. For 
alternative definitions of the d-functions see the selection compiled in [ ] and references therein. Recursion 
formulae and symmetry relations may be used to compute rapidly the Wigner functions in the basis of 
either complex [27,6,17] or real [ , ] spherical harmonics. In the implementations described in this work, 
we employ the recursion formulae described in [ ]. 

The d-functions satisfy a number of symmetry relations. In this work we make use of the symmetry 
relations 

4n(/3) = (-l)'"""rfi„(/3) , (13) 
diJn - P) = i-iy-" dt^JP) (14) 

and 

rfi„(-/3) = (-l)'"-"rfi„(/9), (15) 

which may be found in [33] . In the derivation of our fast spin spherical harmonic transform we often consider 
(i-functions of the form (i^^(7r/2). For this special case, we may infer 

„,„(7r/2) = (-1)^+- „,_„(7r/2) (16) 

from (13) and (14), which implies that (i^^(7r/2) is zero for n = and odd £ -\- m. We make continual use 
of these symmetry relations later in the derivation of our fast algorithms. 

Finally, we note an important Wigner d-function decomposition. The d-function for an arbitrary argument 
may be decomposed into sums of d-functions for a constant argument of 7r/2 [22]: 

i 

dL(/3)=i"-™ E dl,^{n/2)dl,„{7r/2)e^"^'^ , (17) 

which follows from a factoring of rotations. The Fourier series representation of dl^^{p) specified by (17) 
allows one to write the spherical harmonic expansion of a function / on in terms of a Fourier series 
expansion of / on the two-torus T^, with / appropriated extended to this domain (as discussed in more 
detail in section Sec. 3.2). Consequently, (17) is fundamental to the derivation of our fast algorithms. 

2.3. Spin spherical harmonics 

Spin functions on the sphere sf ^ L^(S^, dl](cj)), with integer spin 5 G Z, are defined by their behaviour 
under local rotations. By definition, a spin function transforms as 

sf'iu;) = e-^'^ sfico) (18) 

under a local rotation by where the prime denotes the rotated function. It is important to note that the 
rotation considered here is not a global rotation on the sphere, such as that represented by an element of 
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the rotation group SO (3), but rather a rotation by x the tangent plane at co. The sign convention that 
we adopt here for the argument of the complex exponential in (18) differs to the original definition [21] but 
is identical to the convention used recently in the context of the polarisation of the CMB [39]. 

The spin spherical harmonics syim{^^ ^) form an orthogonal basis for L^(S^, dl](a;)) spin s functions on the 
sphere for \s\ < £. Spin spherical harmonics were first developed by Newman & Penrose [ ] and were soon 
realised to be closely related to the Wigner /^-functions by Goldberg [11]. Spin functions may be defined 
equivalently from the expansion in Wigner /^-functions Df^^{a^ f3^j) of a function in L^(S0(3), dg{p)) 
evaluated at 7 = 0, for fixed n. Consequently, the functions D^^^cp^O^O) or D^_g{(p^O^O), for fixed 5, 
naturally define an orthogonal basis for the expansion of spin s functions on the sphere. The appropriately 
normalised spin spherical harmonics may therefore be given by 



sYe^ie, ^) = (-l)^y Di:_,{^, e, 0) . (19) 

Noting the decomposition given by (11), the spin spherical harmonics may also be defined in terms of the 
reduced Wigner d-functions: 



sYe^ie,^) = {-iy^'^dl^_M e""*^ . (20) 
For completeness, we also state the explicit expression derived by [1 ^ ] for the spin spherical harmonics: 




^ ^ ^ ] (-1)^-^-^ cot^-^+2^(6>/2) , (21) 

where the sum is performed over all values of r such that the arguments of the factorials are non- negative. 
The scalar spherical harmonics are identified with the spin spherical harmonics for spin s = 0, through the 
relation 



rf™oW = ^|:p^^'r(cos^), (22) 

i.e. oYimi^) — Ygrai^uj). The conjugate symmetry relation given for the scalar spherical harmonics generalises 
to the spin spherical harmonics by sYl^{uj) = {—iy~^^ -sY£^_rn{^) • 

The orthogonality and completeness of the spin spherical harmonics follows from the orthogonality and 
completeness of the Wigner D-functions and read 



dQ{ij) sYim{(^) sYi1^,{u) = SwSmm' (23) 

g2 



and 



£ J2 sYU^) ,Y;^{io') = 5{lo - Lo') (24) 

respectively. The spin spherical harmonics form a complete, orthogonal basis for spin functions on the sphere, 
hence any square integral spin s function on the sphere sf ^ L^(S^, drt{u)) may be represented by the spin 
s spherical harmonic expansion 



oo 



£=0 m^-£ 

where the spin spherical harmonic coefficients are given in the usual manner by 

Jem= f dniu;)Jiu;),Y;M- (26) 

JS2 
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The conjugate symmetry relation of the spin spherical harmonic coefficients is given by 
-sfi,-m = (—1)^+^ sfim ^ function satisfying 5/* = -5/ (which, for a spin 5 = function, equates to 
the reality condition) and follows directly from the conjugate symmetry of the spin spherical harmonics. 

Spin raising and lowering operators, S and § respectively, exist so that spin s ± 1 functions may be 
defined from spin s functions [21,11]. When applied to a spin s function the resulting function transforms 
as (dsfYiuj) = e~'^^^~^^'^^{dsf){uj) and {dsf)'{u;) = e~**^^~-^^^(§5/)(cj), where the primed function again 
corresponds to a rotation by x the tangent plane at co. The spin operators are given explicitly by 



and 



sm 



sm 



d id 
dO sin dip 

d id 



sm 



d9 sin 9 dip 



sin' 9 J 



(9, if). 



(27) 



(28) 



3. Fast algorithms 



The harmonic formulation of our fast spherical harmonic transform algorithm for spin functions on the 
sphere is derived here and contains the essence of the fast algorithm. The harmonic formulation involves 
recasting the spherical harmonic transform on the sphere as a Fourier series expansion on the torus. A 
number of subtleties then arise due to the method used to extend a function defined on the sphere to the 
torus and associated symmetry properties. After discussing these subtleties we then present the fast forward 
an inverse spin transforms in detail, including a number of practical issues pertaining to an implementation. 



3.1. Harmonic formulation 

Although we initially wish to compute the spherical harmonic coefficients sfim ^f a spin s function 
sf G L^(S^, drt{uj)), i.e. we wish to evaluate the forward transform specified by (26), we will instead consider 
the inverse transform specified by (25). Using the inverse transform to relate the function 5/ to its spherical 
harmonic coefficients sfim advantage that it doesn't require any quadrature rule to approximate 

an integral over the sphere and, consequently, it is exact, up to numerical precision. It will then be possible 
to construct a fast and exact direct transform by inverting the harmonic representation of the fast inverse 
transform. 

Starting with the representation of the spin spherical harmonics in terms of the reduced Wigner d-functions 
given by (20), and then substituting the d-function decomposition given by (17), the spin spherical harmonics 
can be written in terms of a sum of weighted complex exponentials: 

I ^ 

,r,«(^,<p) = (-l)^r(-+«)y^ ^ dl,^(7r/2)di,,_,(7r/2)e^(-^+'"'^). (29) 

m' ——i 

Substituting this expression for the spin spherical harmonics into the inverse spin spherical harmonic trans- 
form defined by (25), we obtain 

where we assume that s/ is band-limited at L, i.e. sfim — ^^^^^-^^ ^^^^ outer summation may be 
truncated to L. After interchanging the orders of summation, (30) may be written as 

L L 

sf{9,^)= E .^^w e'('"^+'"''') , (31) 

m= — L m' — — L 
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where 



.,F^^, = {-ly M^dl,^{nl2) di,,_,(7r/2) ,f,^ . (32) 

We adopt the convention here and subsequently that ah terms indexed by i and either m or m' are zero for 
|m| > ^ or \m'\ > i. This convention may also be represented explicitly by replacing the lower limit of the i 
summation in (32) from ^ = to ^ = max(|m|, \m'\). 

Notice that (31) is very similar to the Fourier series representation of 5/. However, the Fourier series 
expansion is only defined for functions that are periodic, i.e. for functions defined on the torus. A function 
5/ on the sphere is periodic in (p but not in 6. In order to apply (31) directly it is necessary to extend 
s/, initially defined on the two-sphere S^, to be defined on the two-torus T^. Once s/ is extended in this 
manner an inverse FFT may be used to evaluate sFmm' rapidly from a sampled version of sf • The harmonic 
coefficients sf irn then be computed from sFmm' by inverting (32) to give a fast forward spin spherical 
harmonic transform. Similarly, sf may be reconstructed by computing sFmm' from the harmonic coefficients 
sfirn through (32), followed by the use of a forward FFT to recover a sampled version of sf from sFrnm'- 
This gives a fast inverse spin spherical harmonic transform. By recasting the spherical harmonic expansion 
on the sphere as a Fourier expansion on the torus we may appeal to the exact quadrature of the Fourier 
integral on uniformly sampled grids {i.e. the Shannon sampling theorem). Consequently, our fast transforms 
are theoretically exact (to numerical precision and stability) for a sampling of the sphere compatible with a 
uniform grid on the torus (a corresponding sampling of the sphere is made explicit in Sec. 3.4.1). We address 
the periodic extension of g/ from S^ to T^ next. 

3.2. Extending the sphere to the torus 

When periodicly extending sf from S^ to T^ it is important to ensure that all relations used previously 
are still satisfied on the new domain. Careful attention must therefore be paid to the periodic extension of 
5/, as a naive extension may introduce inconsistency. 

We extend 5/ to sf ^ L^(T^, d^x) by allowing to range over the domain [— tt, tt) and assuming 5/ is 
periodic with period 27r in both (p and (note that d^x = dO dcp is the invariant measure on the (6>, (/:?)-plane). 
Equivalently, we may represent 6> G [0, 27r) and both domains are used interchangeably depending on which 
is most convenient. Representing sf over this extended domain is redundant, covering the sphere exactly 
twice. Nevertheless, such a representation facilitates our fast algorithm, which more than compensates for 
this redundancy. Such an approach is not uncommon, as [9,13] also oversample a function on the sphere in 
the direction in order to develop alternative fast scalar spherical harmonic transforms. 

The method used to extend s/ to the new 27r domain remains to be chosen. Two natural candidates are 
to extend s/ in 6> in either an even or odd manner about zero. We obtain two candidate periodic extensions 
of s/ on T^: s/® and s/°, where the superscripts denote the even and odd extensions respectively, i.e. 
sf^{—0,(p) = sf^{0,(p) and sf'^{—0,(p) = —s/°(6>, (/:?). The appropriate periodic extension to use depends on 
the symmetry properties of the basis functions used to represent the inverse spherical harmonic transform. 
In turns out that we must use both periodic extensions g/® and In Fig. 1 we illustrate these periodic 
extensions of a function on the sphere to give a function defined on the torus. We use an Earth topography 
map to define our initial function on the sphere in this example. 

Let us recall the inverse spin spherical harmonic transform given by (25). To ensure consistency, we should 
extend sf to T^ so that it satisfies the same symmetry properties as sl^m would when extended to this 
domain. For the spin s = case, one might naively expect that Yirn would exhibit even symmetry in 
about zero when extended to 6 e [— tt, tt) based on the definition of Yim given by (1). However, in deriving 
the alternative harmonic formulation of the inverse spherical harmonic transform in Sec. 3.1 we used the 
representation of sYim specified by (20), or equivalently (21). Although all of these definitions are obviously 
identical on the sphere S^, they are not necessarily identical on the torus T^. Since we have used (20) in 
our harmonic formulation we must be sure to satisfy the symmetry properties of this relation on T^. The 
symmetry of (20) in is identical to the symmetry of dl^g{6). Recall from (15) that 
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(a) ,/GL2(S2,d^7M) 



(b) sr e L2(t2, d^a.) 



(c) gL2(t2, d^^r) 



Fig. 1. Periodic extension to the torus of Earth topography data defined originally on the sphere. The left panel shows the 
Earth data defined on the sphere, while the centre panel shows the even periodic extension of the data to the torus and the 
right panel shows the odd periodic extension. Both types of periodic extension are used in our fast algorithm. 



(33) 



i.e. d^g{0) is even for even m -\- s and odd for odd m -\- s. This same symmetry property is also apparent 
from (21). ^ This symmetry could naively appear to be problematic. There is no way we can extend s/ to 
T^ to satisfy this relation. However, we may solve this problem by considering both /® and /°, as alluded 
to earlier. For even m + s, (33) will be satisfied for g/®, whereas for odd m + 5 it will be satisfied for 5/°. 

Our fast spherical harmonic transform then proceeds as follows. We perform two inverse FFTs to compute 
rapidly sFmm' s-^mm' both 5/® and 5/° respectively. We then invert (32) to compute the harmonic 
coefficients slim slim corresponding to 5/® and s/° respectively. The harmonic coefficients sfim of s/ 



will be identical to sfim ^ 

sfim - 



s even only, and the identical to g/^^ for m + 5 odd only, i.e. 




m - 



s even 
s odd 



(34) 



3.3. Symmetry relations 

The Fourier coefficients sFmm' satisfy a number of symmetry relations which may be exploited to increase 
the speed of implementations. Here we examine the symmetry properties of sFmm' using two separate 
approaches. Firstly, we determine the symmetries of sFmm' sFmm' based on the Fourier series represen- 
tations of sf^ and 5/°. Secondly, we examine the symmetries of sFmm' computed directly from sfim using 
(32). Recall that the symmetries of sFmm' will only be satisfied by s^^m' for m + s even and by sFmm' fo^ 
m + s odd. 

Although in practice we will compute sF^^, and sF^m' by taking FFTs of sampled versions of s/® and 
s/° respectively, we may nevertheless represent sF^m' s^mm' by the usual projection of the continuous 
functions on to the complex exponentials, i.e. 

sFmm' = r r ^/(^' ^) e-^(-^+-'^) , (35) 

[ZTry Jo Jo 

where, here and subsequently, we drop the superscript denoting the even or odd function to indicate that 
the relation holds for both (in fact, the discrete Fourier transform merely gives an exact quadrature rule 
for evaluating (35)). By considering all combinations of negative m and we obtain trivially the following 
symmetry relations for ^F^^/from (35): 



^ Here the symmetry of sVgm is dictated by the symmetry of the cot(^/2) term. Since cot is an odd function the overall 
symmetry depends on the power that cot(^/2) is raised to. All cot terms in the sum have either an even or odd power since 
the power increases by 2r for each term. The overall symmetry therefore depends on the parity of s — m, or equivalently s + m. 
Thus sV^m is even for even s + m and odd for odd s + m. 
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-s^-m,m' ~ s^mm' sf — -sf'i (37) 

-.Ff„,=,F^„,* for ,/* = _,/. (38) 
Similarly, we obtain the following symmetry relations for gF^^r. 

sFra^-m' ~ sFmm''') (39) 

-a^V™' = (-l)«i^w* for ./* = _./; (40) 

=,F°^,* for ,/* = _,/. (41) 
Now we derive symmetry relations for sFmm' computed directly from sfim through (32). Firstly, we find 



= (-1)^ ^ ^(7r/2) /_ ^, _,(7r/2) 



= (_1)« r a/^^ (-1)'+™ dl,^i7r/2) (-1)^- rfi,,_.(7r/2) ,/,^ 

— ( 1) sFmm' -) ('^^) 

where in the second line we have applied the d-function symmetry relation (14). For the second symmetry 
relation we find 



£=0 



(r)-(-+«) ^ d^^,^(^/2) d^, _,(7r/2) (-1)™+%/;^ for ,/* = 



47r 

(-l)-+%Fw* for ,/* = _,/, (43) 

where in the second line we have applied the d-function symmetry relations (13) and (14), and noted the 
harmonic conjugate symmetry relation for 5/* = -5/. It is possible to derive the third symmetry relation 
for sFmm' in a similar manner to the previous two, however it is much easier to use these two previously 
derived relations directly, in which case it may be shown trivially that 

-sF-m,-m' = sFmm'* for sf* = -sf- (44) 

Notice that, as expected, we find that the symmetry relations of sFmm' are only satisfied by sF^^, for m + s 
even and by sF^^, for m + 5 odd. 

3.4. Forward and inverse transform 

The harmonic formulation of our fast algorithm was derived in Sec. 3.1. Since this derivation a number 
of practicalities have been discussed, such as extending sf defined on S^ to give sf defined on T^, and the 
symmetry properties of the Fourier coefficients of g/ on T^. We are now in a position to describe our fast 
algorithm in detail. Firstly, we consider the discretisation of s/ on and the inversion of (31) using an 
FFT to compute sFmm' rapidly. We then consider the inversion of (32) to compute the spherical harmonic 
coefficients sfirn- Finally, we discuss how our formulation may also be used to compute the inverse spherical 
harmonic transform, i.e. to recover 5/ from the harmonic coefficients sfim- Although our fast forward 
and inverse spin spherical harmonic transforms appear to follow directly from (31) and (32), a number of 
subtleties arise due to the symmetries of the functions involved. 

3.4.1. Computing Fourier coefficients on the torus 

We first discretise sf defined on the two-torus T^. For now we are not concerned with the method by 
which s/ is periodicly extended to T^, i.e. through either even or odd symmetry, hence we drop the related 
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superscript. Previously, we extended 5/ to by extending the 6 domain to [— tt, tt). However, we noted 
that sf could equivalently be represented on the domain 6> G [0, 27r). For convenience, we choose this second 
domain here since it simplifies the subsequent discretisation. We adopt the equi-angular sampled grid with 
nodes give at 

^t = ^^|^, wheret = 0,l,...,2L (45) 

and 

^P = ^^^, wherep = 0,l,...,2L. (46) 

An odd number of sample points are required in both 6 and (p so that a direct association may be made 
between the sampled 5/ and sFmm' (where both m and m' range from — L, . . . , 0, . . . , L). Moreover, we make 
the association = 2L + 1, where N is the number of samples in both the 6 and Lp dimensions, to ensure 
that Nyquist sampling is satisfied. The node positions specified by (45) and (46) eliminate repeated samples 
at the poles ^ = and = 27r, since these points are excluded from the grid. However, it is not possible to 
eliminated repeated samples at ^ = tt, since we require a discretisation that is symmetric about tt but which 
contains an odd number of sample points. The sampled version of s/ is given by 

shpA=shOu^p). (47) 

where we use the convention that square brackets are used to index the discretised function. 

Now that we have sampled s/, we may invert (31) using an FFT. This process is simple, but not trivial, 
and we discuss the various approaches that may be used to solve this problem. We start by discretising (31): 

s~f\pA= E E .F«™,e-M2p+i)W(2*+i)]/(2L+i)_ (48) 

m——L m' ——L 

It is necessary to rewrite (48) in a form that is suitable for the application of standard FFT algorithms, 
where we adopt the forward discretise Fourier transform definition 

N-l N-1 

Gw = H9\P, ^]} = E E 5b' *] e-*2-('"P+'"'*)/^ , (49) 
with inverse 

N-l N-l 

5b, t] = T-'{G^^,} = ^ E E e'2-(-f+-'*)/^ , (50) 

m=0 m'=0 

for indices p, t, m, m' all ranging over 0, 1, . . . , A" — 1. In order to convert (48) into a form similar to (50) 
so that an FFT may be applied to evaluate it, we first shift the summation indices to give 

2L 2L 
m—0 m'—Q 

from which it follows that 

gi2.L(p+*+l)/(2L + l) J^^^^ ^ (2L+ 1)2^-1 {e^^('"+'"')/(2i + l) sFra-L,ra'-L} ■ (52) 

This may be easily inverted using the definition of the discrete Fourier transform, giving 

sFm-L,m'-L = ^ ^|,i2.L(p+*+l)/(2L + l) J^^^^^ (53) 

This may be evaluated through an FFT directly, applying the appropriate phase shifts at each stage to recover 
sFm-L,m'-L- Alternatively, we may represent the phase shift applied to sf[p^t] in the spatial domain as a 
spatial shift in the frequency domain. In this case we recover 

-i7r(m+m')/(2L+l) . . 
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where m and now range from — L,...,0,...,L and the FFT-shift operator S interchanges quadrants 
of a sampled two-dimensional signal by interchanging quadrant I with IV and also interchanging quadrant 
II with III. The final representation given by (54) is marginally less computationally demanding than the 
representation given by (53) since it is not necessary to compute and apply one of the phase shifts. Moreover, 
(54) involves an FFT of a real signal for real spin s = functions, rather than a complex one, for which 
FFTs are approximately twice as fast and require half the memory requirements. In our implementation we 
use the expression given by (54) to compute sFmm' rapidly using an FFT. Finally, we note that (54) may 
be inverted easily through 

J\p,t] = (2L+ 1)2 ^-i5-i|e-(™+™')/(2L+i) j^^j ^ (55) 

where is the inverse FFT-shift operator. 

We have so far ignored the fact that sf always exhibits either even or odd symmetry. This property may 
be exploited by using fast cosine or sine transforms in place of the FFT, which will improve the speed of the 
algorithms but will not alter their complexity. For now, we ignore this optimisation but intend to readdress 
it in future. 



3.4.2. Forward transform 

Now that we have discussed how to compute sFmm' rapidly from a sampled version of 5/, we turn our 
attention to inverting (32) in order to recover the harmonic coefficients sfim- Recall that (32) is not satisfied 
simultaneously by the two periodic extensions of sf on T^, i.e. sf^ and for all m -\- s. Instead, it is 
satisfied by sf^ for even m-\- s and by s/° for odd m-\- s. To recover the harmonic coefficients sfirn^ we solve 
(32) for sfim even m -\- s and solve for sfim for odd m -\- s. 

Let us consider fixed m = mo but all values of m\ Throughout we also consider fixed s only. For this case 
we may rewrite (32) as the matrix equation 

/ /2L + 1 \ 



/ sJ^ rao^ — L \ 



s-^mo ,0 



\ sFmo,L I 



-(mo+s) 



47r 



4^ "0,mo "0,-s 



^2L + 1 


^-l,mo 




47r 


/2L + 






V 47r 


/2L + 


^d^ 




V 47r 



v 



2L 



47r 



( -fo,mr^ 



\sfL,mo J 



J 



where we have used the convention that d-functions with no argument specified take the assumed argument 
of 7r/2, i.e. df^^ = dl^g{7r/2)^ and empty elements of the matrix not enclosed by ellipsis dots are zero. In 
order to solve this equation for the harmonic coefficients fimo ? we only require the lower half of the system; 
the upper half of the system will be satisfied automatically by symmetry considerations. We thus consider 
only the triangular system 



mo ,0 



-(mo+s) 



\sFmo,L I 



d^ d^ 

4^ «0,mo ^0, 




■d'^ d'^ 
4^ "0,mo "0,- 



^-^d^ d^ ^ 

4^ "'0,mo "'0,-s 



,4- 



2L + 1 
47r 



2L 



47r 



d^ d^ 

"l,mo "1,- 



^L,mo ^L,-s / 



s/0,mo 



\sfL,mo J 
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which we write 



s'^mo sf n 



(56) 



where the ^-(^o+s) factor is incorporated into the definition of s^mo- Once the sFmm' coefficients are 
computed, as outlined in Sec. 3.4.1, we can then compute sfimo for all £ by inverting the triangular system 
specified by (56). For an arbitrary complex signal 5/, a similar system must be solved for all mo- However, 
if one is computing for both spin ±s or for 5 = for a function sf that satisfies g/* = -sf {i.e. a real signal 
for the spin 5 = case), then it is only necessary to solve systems for non-negative mo. Conjugate symmetry 
may then be used to compute the harmonic coefficients for negative mo- 

For the spin s = case, additional symmetries exist in the matrix sD^^ that may be exploited to improve 
the efficiency of inversion. Recall that dl^Q{7T/2) is zero for odd i-\-m (which follows from (16)). This implies 
that all entries of oDmo ^ire zero for odd £ + m^ We illustrate this point with an example for L = 5. In this 
case the matrix oDmo takes the form 



• • • 



oDmo 



V 



where all entries without a dot are zero. Due to the structure of oDmo it is possible to separate (56) into 
two independent triangular systems by extracting out the even and odd rows of the original system. This 
provides a small improvement in the number of computations required to invert (56) and is the method 
employed in our implementations for the spin s = case. 

We have so far discussed how to invert (32) for sf on T^ in order to recovered sfem- briefly 
review how to recover the harmonic coefficients sfim of the original signal sf on S^. Recall that (32), and 
consequently (56), are satisfied by g/® and 5/° only for even and odd m-\- s respectively. Consequently, the 
harmonic coefficients of sf will be given by 



sfim — 




m-\- s even 
m -\- s odd 



(57) 



3.4.3. Inverse transform 

Our fast algorithm may also be used to perform the inverse spherical harmonic transform, i.e. to recover 
a function on the sphere from its spherical harmonic coefficients sfirn- ^^^^ ^^^^ cannot construct 
sfim sfim from sfim siucc wc do uot kuow what values sfim sfim should take for odd and even 
m -\- s respectively. However, this is not problematic since we can reconstruct an 5/ that is extended to T^ 
in a different manner, but which is identical to 5/ only on the half of T^ that maps directly onto S^. It is 
possible to do this by ffist applying (56) on the harmonic coefficients sfim directly. Values for negative m' 
may then be computed using the symmetry relation specified by (42). An FFT may then be used to recover 
a sampled version of sf through (55). We then recover sf from the values of sf ^oi which G [0,7r). Note 
that outside of this domain sf may not have even or odd symmetry, however we are not interested in these 
values and thus discard them. 
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4. Complexity and numerical experiments 

The fast forward and inverse spin spherical harmonic transforms derived in the previous section have been 
implemented using the FFTW^ library to compute Fourier transforms. We discuss here the theoretical 
memory and computational requirements of these algorithms, before performing numerical experiments to 
evaluate execution time and numerical stability To avoid unnecessary complexity we restrict the numerical 
experiments to real spin s = functions on the sphere. 

4.1. Memory and computational complexity 

It is possible to reduce the computational burden of our fast transforms by precomputing the reduced 
Wigner (i-functions for an argument of 7r/2 over a range of ^, m and n values. This type of precompute 
compares favourably with the precompute required for other fast scalar spherical harmonic transforms, 
which require a precomputation of associated Legendre functions for a range of i and m indices, and over 
the values of a discretised grid. Our precompute of the d-functions is independent of the grid: as the band- 
limit L increases we must precompute more terms, but it is not necessary to re-evaluate the precomputed 
terms already calculated for a lower band-limit. This is not the case when precomputing associated Legendre 
functions over the 6 grid. When evaluating the d-functions it is only necessary to precompute df^^{7r/2) for 
no n- negative m and n. Values for negative m and n indices may then be recovered by noting the symmetry 
relations (13) and (16). Naively terms must be computed, however this ignores the restriction that \m\ < i 
and \n\ < ^ for each i. Taking this restriction into consideration the number of precomputed terms is 

L 

^ (£ + 1)2 = + 3L^/2 + 13L/6 + 1 . 

To precompute all of the associated Legendre terms (for m > only and noting the restriction that \m\ < £ 
for each £) requires the computation of terms, when is defined over a grid of size 2L. Although the 
complexity of both precomputes scale as (9(L^), precomputing the reduced Wigner d-functions requires 
approximately one third of the number of terms as precomputing the associated Legendre functions. 

The overall memory requirements of our fast spin spherical harmonic transform algorithms are dominated 
by the storage of the precomputed reduced Wigner d-functions and therefore scale as 0{L^). This precom- 
puted data could be read from disk rather than stored in memory, however this would reduce the speed of 
the algorithms considerably. In this setting the memory requirements scale as 0{L'^). 

We now consider the computational complexity of our algorithms. Although we have discussed the details 
of the algorithms in Sec. 3.4, the algorithms essentially equate to directly evaluating or inverting (31) and (32) 
as derived in Sec. 3.1. All other operations do not alter the overall complexity of the algorithms and merely 
introduce a prefactor. By applying an FFT it is possible to reduce the complexity of evaluating/inverting 
(31) from 0{L^) to 0{L'^ log2 L). The linear system defined by (32) involves a triangular matrix, as specified 
explicitly by (56). A triangular system of size L x L may be evaluated/inverted in 0{L'^) operations. This 
system must be considered for each value of mo, hence the overall complexity of evaluating or inverting (32) 
is 0{L^). The complexity of both of the forward and inverse fast spin spherical harmonic transforms is thus 

^ V ^ V ^ V " 

Evaluation/inversion of (31) by FFT Evaluation/inversion of (32) Overall complexity 

Using our implementation of these algorithms with precomputed data stored in memory, we evaluate 
computation times for real spin s = signals over a range of band-limits. These timing tests are performed 
on random band-limited test functions on the sphere. The test functions are defined through their spherical 
harmonic coefficients, with independent real and imaginary parts distributed uniformly over the interval 
[— 1, 1]. Six band-limits from L = 16 to L = 512 are considered, increasing by a factor of two at each step. 
Computation time measurements are performed on a laptop with a 2.2GHz Intel Core Duo processor and 

^ http://www.fftw.org/ 
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Fig. 2. Average computation time r for the forward and inverse transforms of random real spin s = test functions. Note that 
the computation time scales like as expected. 

2GB of RAM and are averaged over five random test functions. Fig. 2 siiows a log-log plot of the average 
computation time r with band-limit. Also shown on Fig. 2 is the slope corresponding to computation time 
that scales as L^. The average computation time for the implementation of our fast spin spherical harmonic 
transforms scales approximately as L^, as predicted by the theoretical computational complexity discussed 
above. 



4.2. Numerical stability 

The numerical stability of our algorithms is examined in this section. Unfortunately the algorithms are 
found to be stable only for very low band-limits in the order of L :^ 32 and below. We first evaluate the 
stability of the algorithms in detail, before examining the source of the instability and possible solutions. 

To test the stability of our algorithms we perform the following numerical tests, restricted here to real 
spin s = functions on the sphere. Firstly, we consider a random test function on the sphere defined 
through its harmonic coefficients, with independent real and imaginary parts distributed uniformly over the 
interval [—1,1] (the same test functions considered in Sec. 4.1). By defining the test function explicitly in 
harmonic space we ensure that it is band-limited. We then perform the inverse spherical harmonic transform 
to compute the associated real values of the function on the sphere, before performing the forward spherical 
harmonic transform to reconstruct the spherical harmonic coefficients of the function. The error between 
the original spherical harmonic coefficients f^rn and the reconstructed spherical harmonic coefficients g^rn 
is examined to assess the numerical stability of the algorithms. Five different error metrics are computed, 
including the maximum absolute error 

Cmax max Ifim - 9lm\ , (58) 

the mean absolute error 

emean = mean \firn - gim \ , (59) 

the median absolute error 

^median = median l/^^ - gim\ , (60) 

the root-mean-square error 

Crms = ^ \fim " gim\ (61) 
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(a) All band-limits (b) Lower band-limits 



Fig. 3. Average errors between the original and reconstructed spherical harmonic coefficients of random real spin 5 = functions 
for the numerical stability experiment explained in the text. 

and the relative root-mean-square error 

erelrms = ^ \fim - gimf / ^ \fimf • (62) 

Error measurements for seven band-limits from L = 8toL = 512, increasing by a factor of two at each step, 
averaged over five random test functions on the sphere are illustrated in Fig. 3. For very low band-limits the 
errors are on the order of the machine precision, however as the band-limit increases the instability of the 
algorithms become apparent as the errors quickly increase to unacceptable levels. From the median error it 
is apparent that for band-limits below L :^ 64 many of the spherical harmonic coefficients are well recovered 
and only a relatively small number of coefficients are responsible for the large errors measured by the other 
error metrics. However, above this band-limit even the median error increases to an unacceptable level 
indicating that most coefficients are poorly recovered. To examine the location of the harmonic coefficients 
that are poorly recovered we plot in Fig. 4 the original harmonic coefficients fim and the reconstructed 
harmonic coefficients girn ^oi one test function for various band-limits. For L = 16 and L = 32 (and L = 8 
although not shown) all harmonic coefficients are recovered accurately. For L = 64 many coefficients are 
recovered accurately, however a small number of coefficients near the diagonal are not. Most coefficients are 
recovered poorly for band-limits above L :^ 128. Errors are first introduced along the diagonal where i ^ m 
and then quickly spread as i increases. 

The source of the numerical stability seen above has been traced to the system defined by (32), which is 
specified explicitly as a matrix system in (56). As an example, in Fig. 5 we show the triangular matrix oDmo 
defining this system for L = 50 and mo = 38. The matrix is extremely poorly conditioned with condition 
number 6 x 10^^. We have investigated a number of numerical methods for solving these ill-conditioned 
systems, including a singular value decomposition (SVD), preconditioning, the generalised minimal residual 
method (GMRES), iterative refinements and various re- normalisations; however, without success. It may 
indeed be possible to solve these systems with a suitable numerical method that we have not yet tried and 
so we make the system shown in Fig. 5, with correct o-^mo and of mo vectors, available for download^ in 
case others wish to examine the system and have more success in solving it. Alternatively, it may be possible 
to rewrite the original formulation of our fast algorithms so that the equations are normalised in such a way 
as to avoid these ill-conditioned matrices, however we have as yet been unable to do this. Nevertheless, the 
conditioning of these systems is only a problem when inverting the system, thus it is only the fast forward 



^ http : www . mrao . cam . ac . uk/"^ j dm57/ research/f sht . html 
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(a) Original L = 16 



(b) Reconstructed L = 16 




20 25 30 



(e) Original L = 32 



(f) Reconstructed L = 32 




(i) Original L = 64 



(j) Reconstructed L = 64 



(c) Original L = 128 



^1 



(d) Reconstructed L = 128 




100 150 200 250 



(g) Original L = 256 (h) Reconstructed L = 256 




(k) Original L = 512 



(£) Reconstructed L = 512 



Fig. 4. Absolute value of original and reconstructed spherical harmonic coefficients for the numerical experiments explained in 
the text. The colour scale ranges over the allowable values of the original harmonic coefficients, i.e. over [0,v^]. 

spin spherical harmonic transform that suffers from this instabihty problem. The inverse spin spherical 
harmonic transform is likely to be accurate to high precision up to band-limits where the reduced Wigner 
(i-functions can be computed accurately. However, we have not tested the inverse transform in isolation since 
a stable forward transform that is exact (to numerical precision) does not exist currently on the (^^, (pp) grid 
positions defined by (45) and (46). 



5. Conclusions 



Algorithms have been derived to perform the forward and inverse spherical harmonic transform of functions 
on the sphere with arbitrary spin number. These algorithms involve recasting the spin transform on the 
sphere as a Fourier transform on the torus. FFTs are then used to compute Fourier coefficients, which are 
related to spherical harmonic coefficients through a linear transform. By recasting the problem as a Fourier 
transform on the torus we appeal to the usual Shannon sampling theorem to develop spherical harmonic 
transforms that are theoretically exact for band-limited functions, thereby providing an alternative sampling 
theorem on the sphere. 
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Fig. 5. Absolute values of the ill-conditioned triangular matrix of the system that must be solved for L = 50 and mo = 38. The 
matrix is extremely poorly conditioned with condition number 6 x 10^^. 

Previous approaches to performing spin transforms on the sphere have appealed to the spin lowering and 
raising operators to relate spin transforms to scalar transforms. For spin s transforms based on this method 
the lowest computational complexity that can be attained scales as 0{sL'^ \og2 L)^ whereas our algorithm 
scales as 0{L^) for arbitrary spin number s. Moreover, our algorithm is general in the sense that the spin 
number is merely a parameter of the algorithm; altering the spin number does not alter the structure of the 
algorithm. Furthermore, the algorithms apply for arbitrary band-limit L, whereas many other methods are 
restricted to band-limits that are a power of two. 

Our algorithms have been implemented. The computational speed of the implementation may be improved 
in future by using fast discrete cosine and sine transforms in place of FFTs, although this will not alter the 
scaling of computation time with band-limit. Numerical tests showed that computation time scales as 
as expected. Numerical tests were also performed to access the stability of the algorithms. Unfortunately 
the algorithms were found to be unstable for band-limits above L 32. The source of the instability 
was determined to be due to the poor conditioning of the linear system relating the Fourier and spherical 
harmonic coefficients. Consequently, we expect only the forward transform to be unstable and not the inverse 
transform, however it is not possible to verify this hypothesis since a stable forward transform that is exact 
(to numerical precision) does not exist currently on our pixelisation of the sphere. A number of numerical 
methods have been applied in an attempt to solve the ill-conditioned system but these have not proved 
successful. 

A fast algorithm to perform a scalar spherical harmonic transform has been developed by Anthony Lasenby 
(private communication) in the context of geometric algebra. This algorithm shares many of the properties 
of the algorithms derived in this work (including, unfortunately, hints of numerical stability issues) and 
may simply be an alternative formulation, for the spin s = case, of the algorithms derived herein. If 
this indeed proves to be the case, then this alternative geometric algebra formulation may provide insight 
on alternative approaches to reformulate our algorithms to improve numerical stability. In ongoing work 
we are investigating the numerical stability problem in more detail and are attempting to renormalise or 
reformulate the algorithms in such a way as to eliminate this instability. 
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